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O : ABSTRACT 

. The thermonuclear runaway that culminates in the explosion of a Chandrasekhar mass white 

dwarf as a Type la supernova begins centuries before the star actually explodes. Here, using a 
3D anelastic code, we examine numerically the convective flow during the last minute of that 
runaway, a time that is crucial in determining just where and how often the supernova ignites. 
c/3 , We find that the overall convective flow is dipolar, with the higher temperature fluctuations in 

an outbound flow preferentially on one side of the star. Taken at face value, this suggests an 
asymmetric ignition that may well persist in the geometry of the final explosion. However, we 
also find that even a moderate amount of rotation tends to fracture this dipole flow, making 
ignition over a broader region more likely. Though our calculations lack the resolution to study 
the flow at astrophysically relevant Rayleigh numbers, we also speculate that the observed dipolar 
flow will become less organized as the viscosity becomes very small. Motion within the dipole 
flow shows evidence of turbulence, suggesting that only geometrically large fluctuations (~ 1 
km) will persist to ignite the runaway. We also examine the probability density function for the 
temperature fluctuations, finding evidence for a Gaussian, rather than exponential distribution, 
which suggests that ignition sparks may be strongly spatially clustered. 

Subject headings: Supernovae, hydrodynamics 

1. INTRODUCTION cused upon the computationally challenging prob- 
lem of how a nuclear fusion flame, once formed 
The leading model for a Type la supernova is in the white dwarf mter ior, succeeds in burning 
a carbon-oxygen white dwarf that accretes mat- a sufficient fraction of its mass to 56 Ni and in- 
ter from a binary companion at a sufficiently high termediate mass elements to power a healthy ex- 
rate that it is able to grow to the Chandrasekhar plosion and giye a credible light curve and spec _ 
mass and explode (eg., Hillebrandt & Niemeyer tmm ( eg ^ Gamez0 et aL 2 003; Plewa et al. 2004; 
2000). Much of the research in recent years has fo- R5pke & Hillebrandt 2005). Equally important, 
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however, is just how the flame ignites. A flame 
starting from a single point at the middle, a single 
point off center, and a multitude of points isotrop- 
ically distributed around the center will all give 
qualitatively different results (Nicmcycr & Hillc- 
brandt 1995; Niemeyer, Hillebrandt, & Woosley 
1996; Plcwa ct al. 2004; Garcfa-Senz & Bravo 2005; 
R6pke et al. 2005). The first numerical simula- 
tion of carbon ignition suggested central ignition 
(Hoflich & Stein 2001), but may not have correctly 
represented the pre-explosive convective flow be- 
cause of its restricted geometry. Analytic calcu- 
lations have supported the idea of off-center igni- 
tion, probably at multiple points (Garcia-Senz & 
Woosley 1995; Woosley, Wunsch, & Kuhlcn 2003; 
Wunsch & Woosley 2004), but require their own 
assumptions, for example the persistence of fluc- 
tuations and the existence (or non existence) of an 
ordered background flow. 

The numerical simulation of carbon ignition 
in this environment involves a different sort of 
physics, and optimally, a different sort of computer 
code than the flame propagation problem. The 
fluid motion is very subsonic (Mach number less 
than about 0.01) and the density contrasts that 
must be followed are small ~ 10~ 5 . Both can pose 
problems for compressible hydrodynamics codes. 
The simulation must be followed in 3D for a long 
time and this is challenging for a code that is 
Courant limited by the speed of sound. However, 
these same circumstances are very favorable for 
anclastic hydrodynamics, which is Courant lim- 
ited by the fluid velocity. In § 3.1 we describe the 
implementation for the supernova problem of a 3D 
anclastic code that has previously been employed 
to study the sun (Glatzmaier 1984), the earth's 
geodynamo (Glatzmaier & Roberts 1995), and 
planetary convection (Sun, Schubert, & Glatz- 
maier 1993). This code is based on spectral meth- 
ods that greatly improve the effective resolution 
over ordinary spatial grids. As we shall see, res- 
olution is a big issue since the Reynolds number 
in the star is Re <~ 10 14 (Woosley, Wunsch, & 
Kuhlcn 2003), far beyond what can be achieved in 
any presently existing code. We are thus only able 
to glimpse large scale features (§ 4.2) and hints of 
complex structures beneath. 

Nevertheless, we come to some interesting con- 
clusions. The overall convective flow, for the range 
of Reynolds numbers that is accessible, is dipo- 



lar. This has been seen before in different environ- 
ments (eg., Woodward et al. 2003), but the present 
study is the first to find it in numerical simulations 
of Type la supernova ignition. This flow suggests 
the supernova may ignite in a lop-sided fashion, 
but there are many caveats (§ 5). We are also 
able to determine the probability density function 
(PDF) of the temperature fluctuations and show 
that it is Gaussian. This has important implica- 
tions for the number of ignition points (Woosley, 
Wunsch, & Kuhlen 2003). We also call attention 
to the turbulence that exists in the outgoing "jet" 
of the dipole (§ 4.4). Superimposed on the flow 
which carries the sparks that will ignite the su- 
pernova is the beginning of what looks to be a 
Kolmogorov cascade to smaller scales. This will 
limit the size of the temperature fluctuation that 
is necessary to ignite the star off center. Sparks 
that are too small will be dissipated by turbulence 
before they run away. It will take a finer resolved 
study than the present one to conclusively address 
the issue, but it may be that the smallest ignition 
points are of resolvable size, ~ 1 km. 

2. The Initial Model 

In our implementation of the anelastic approxi- 
mation we solve for thermodynamic perturbations 
to an iscntropic, constant reference state (see Sec- 
tion 3.2). The one-dimensional background state 
has been determined here using an implicit hydro- 
dynamics code, kepler( Weaver, Zimmerman, & 
Woosley 1978). A white dwarf composed of 50% 
carbon and 50% oxygen was allowed to accrete 
(carbon and oxygen) at a rate of 10~ 7 M Q yr _1 
When the central density and temperature 
reached 3.2 x 10 9 g cm" 3 and 2.5 x 10 8 K, nuclear 
energy generation from carbon fusion exceeded the 
plasma neutrino loss rate. The excess energy was 
carried away by convection. At this point the 
star's mass was 1.38 M© and its radius 1580 km. 

Over the next 5000 years the central temper- 
ature rose and the extent of the convective core 
grew, from initially less than 1%, to 0.45 M when 
the central temperature was 4 x 10 8 K (2 years be- 
fore explosion), and to 0.95 M when the central 
temperature was 6 x 10 8 K (3 hours before explo- 
sion). During this ramp up to runaway, neutrino 
losses and radiation transport to the surface of the 
star were negligible and a small amount of energy 
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Fig. 1. — The one-dimensional kepler model: Temperature, density, pressure, entropy per baryon, nuclear 
energy generation, and radiative luminosity as a function of radius. Thick lines are for the T c ,8 = 7, thin 
lines for the T c $ — 7.5 snapshot. The vertical dotted lines denote the inner and outer boundary of the 
three-dimensional anelastic simulation. Dashed lines indicate the iscntropic reference state (see Section 3.2). 
The gray region in the entropy plot demarcates the range of radii that are convectively unstable in the 
kepler model. 



went into expansion (the central density first in- 
creased slightly then declined by 15%). Most of 
the energy went simply into raising the tempera- 
ture of the convective core and extending its extent 
(Baraffe et al. 2004). 

When the central temperature reached 7.5 x 10 8 
K, the convective core was 1.15 M Q (1040 km) and 
the central density, 2.7 x 10 9 g cm~ 3 . The net 
binding energy of the white dwarf (internal plus 
gravitational) was 4.36 x 10 50 erg and the pressure 
scale height was 450 km. If mixing length convec- 
tion was left on, the explosion occurred 190 s later. 
Radial profiles of T, p, P, sm p /kB (entropy per 
baryon), e nuc , and L are shown in Figure 1, at two 
different times: when the central temperature has 
reached T c .$ = 7 and T c .$ = 7.5. It is interesting 



to note that though the entropy in the convective 
region in the Kepler model was constant to within 
1%, the convective flux was not constant for any 
appreciable span of radii. Two-thirds of the total 
nuclear energy (3.0 x 10 45 erg s _1 ) was produced 
within the inner 0.01 M Q , but the luminosity it- 
self rose to a maximum of 2.6 x 10 45 erg s _1 and 
declined to approximately zero at the outer edge 
of the convective core. Almost all of the nuclear 
energy is lost in transit to heating along the way. 
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3. A THREE DIMENSIONAL MODEL 
USING ANELASTIC HYDRODYNAM- 
ICS 

3.1. The Code and its Modifications 

The code used for the simulations presented 
here is based upon the anelastic hydrodynamics 
code of Glatzmaicr (1984). The hydrodynamic 
equations are implemented in a fully spectral man- 
ner, expanding all quantities in spherical harmon- 
ics to cover the angular variations, and in Cheby- 
shev polynomials for the radial dependence. Al- 
though a spectral method in spherical coordinates 
is more difficult to implement than one based on 
a finite differencing scheme on a Cartesian grid, 
and also requires more communication on a multi- 
processor parallel computation system, a major 
advantage is the greater efficiency provided by 
the spectral method. To achieve a modest ac- 
curacy of a few percent, spectral methods typi- 
cally require only half as many degrees of free- 
dom per dimension fourth order finite dif- 
ference method (Boyd 2000). Chebyshev poly- 
nomials are defined on a domain from —1 to +1, 
which in this code is mapped to the inner and 
outer boundary of the convective region. Unfor- 
tunately, the current implementation of spherical 
coordinates prohibits the computational domain 
from extending all the way to the origin, since 
1/r terms in the angular derivatives lead to di- 
vergences. A central sphere is cut out and re- 
placed by an inner impermeable boundary condi- 
tion, even in situations where the convective region 
includes the center of the star (see 3.3). The code 
employs a rhomboidal truncation of the spherical 
harmonic modes, which means that every longi- 
tudinal Fourier mode (e tm ^) has the same num- 
ber of latitudinal associated Legendre functions 
(P z m (cos0)). A rhomboidal truncation scheme is 
easier to parallelize than a more conventional tri- 
angular truncation (—1 < m < +1), and it ensures 
that the latitudinal resolution is the same for all 
longitudinal modes. 

In this formulation of anelastic hydrodynamics, 
the mass, momentum, and entropy conservation 
equations are expanded in a power series around 
a one-dimensional, constant, and isentropic refer- 
ence state, and only the lowest orders terms are 
retained. The reference state is determined by 
the one-dimensional stellar evolution code kepler 



(see Sections 2 and 3.2). In the following all refer- 
ence state quantities are barred, and all perturba- 
tions denoted with a 5. The fundamental quanti- 
ties in this code are the entropy (5s) and pressure 
(5P) perturbations. In order to relate these to the 
more familiar quantities density and temperature, 
we require an equation of state. 



ST 



dT 



5s 



dT 
dP 



5P (la) 



(ib) 



The equation of state is specified by four ther- 
modynamic partial derivatives, which are part of 
the time-independent reference state and depend 
only on radius. We use the Hclmholtz equation 
of state code (Timmes & Swesty 2000) to calcu- 
late these derivatives for a given KEPLER reference 
state. 

In total, we solve for two thermodynamic vari- 
ables (5s and 5P), the perturbation of the gravita- 
tional potential (5^ g ), and the three components 
of the fluid flow velocity (v r , vg, and v^). The rel- 
evant equations are: 



V • (pv) 



(2a) 



,5P 



-V • (pv v) - /oV(— + 5<P g ) 



dp 



5sgr + 2pvxft 



+V-(2pP(E--(V-v)l)) (2b) 



-85s 
~dt 



V • (RpTVSs) - TV • (p5sv) 
+p(e nuc + 5e 

-p(e v + 5e v ) (2c) 
V 2 5$ g = AnG6p (2d) 



Here E is the rate of strain tensor, I the identity 
matrix, and f2 the reference state rotation vector. 
In order to allow for the differential heating of hot 
and cold convective elements, we have included 
a first order perturbation term in addition to the 
one-dimensional background nuclear burning rate. 
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This background energy generation rate is deter- 
mined from the KEPLER output and fitted to a 
power law, pT a . For the simulations discussed in 
this paper a is approximately equal to 27. The 
total energy generation rate is then: 



Sp ST s 

< nac = Enuc + ^uc « P T ° ( 1 + ~= + a Tj7 ) ^ 



Energy losses to neutrinos {e v , e v ) act as a sink 
in the entropy conservation equation. Results 
from one dimensional simulations indicate that 
these neutrino losses are negligible (Section 2), and 
we have included them in Equation (2) only for 
completeness. 

Momentum and energy transport at scales be- 
low our spatial resolution are handled by introduc- 
ing artificially high viscous and thermal diffusion 
coefficients v and R. In our simulation the Prandtl 
number (Pr = D/R) is kept at unity, and v and R 
are lowered as much as possible for a given reso- 
lution (see Section 3.3). These turbulent diffusion 
coefficients also enter in the calculations of the di- 
mcnsionless numbers characterizing the fluid flow, 
the Rayleigh (Ra) and Reynolds (Re) numbers. 



Ra 
Re 



D 3 gA& 
CpDR 
vD 



(4a) 
(4b) 



Here D is the depth of the modeled region and 
As the change in specific entropy across D. Ra 
is proportional to the ratio of buoyancy to dif- 
fusion forces, and Re is proportional to the ra- 
tio of inertial to viscous forces. A large Rayleigh 
number indicates very vigorous convection, and a 
large Reynolds number indicates a highly turbu- 
lent flow. Typical Rayleigh numbers in the cen- 
tral convective region of a pre-SNIa white dwarf 
are around 10 25 and typical Reynolds numbers 
around 10 14 (Woosley, Wunsch, & Kuhlen 2003). 
Use of the turbulent viscous and thermal diffu- 
sion coefficients instead of the true molecular ones 
greatly reduces the magnitude of Ra and Re that 
are achievable with this code. Our simulations 
reach Ra s=s 10 7 and Re w 1500, just entering the 
turbulent convective regime, but a far cry from the 
true physical conditions (see Section 4). 



3.2. Mapping the Background State into 
the 3D-Code 

From a given kepler output (Section 2) we 
extract the density and pressure profiles for the 
region of interest. For numerical stability, we then 
fit a polytrope (P(r) = Kp(r) 1+1 ' n ) to this pro- 
file, and solve the resulting Lane-Emden equation 
for P(r) and p(r). Using the Helmholtz equa- 
tion of state code (Timmes & Swesty 2000) we 
solve for the corresponding temperature profile 
T(r) and for the thermodynamic partial deriva- 
tives (see Section 3.1). Since the kepler models 
are close to, but not completely isentropic we also 
determine a volume averaged background state en- 
tropy. This quantity is useful for comparisons with 
analytical models, but never used in the 3D anelas- 
tic calculation itself. 

3.3. Models and Procedures 

Analytical considerations show (Woosley, Wun- 
sch, & Kuhlen 2003) that the ignition of the fi- 
nal deflagration leading to the supernova explo- 
sion is likely to start somewhere off-center at about 
150 — 200 km. We have chosen to model a region 
extending from 50 km to 500 km. The bound- 
ary conditions at these radii are impermeable and 
stress-free (v± = and f (^) = 0). 

Since almost all of the energy generation and 
the most vigorous convection occur at radii less 
than 200km , we don't expect this outer bound- 
ary at 500km to unduly influence our results. The 
central boundary is more problematic, since in a 
real star nothing prevents fluid from streaming 
through the center. In our model cold sinking fluid 
will splash around the inner boundary. Only by 
extrapolating from conditions outside R = 50km , 
can we say anything about the flow pattern right 
at the center. 

Additional boundary conditions are necessary 
for the entropy equation. Since the radial veloc- 
ity, and with it, the convective heat flux, is forced 
to zero at the inner and outer boundary, turbu- 
lent diffusion is the only mechanism by which en- 
ergy can be transported across these boundaries. 
This heat diffusion is carried by a negative ra- 
dial entropy gradient. At the inner boundary the 
kepler luminosity is essentially zero, and so we 
set dSs/dr = there. At the outer boundary 
the KEPLER luminosity is non-zero (~ 10 45 - 10 46 
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Table 1: Summary of models 



Name 


-*c,8 ^max 


^ max I m ax A t 


n 


Description 






Is] 


[rad s- 1 ] 




C7 


7.0 241 


42 85 70.36 


0.0 


fiducial model 


C7rot 


7.0 


77.06 


0.167 


rotating background state 


C75 


7.5 


40.90 


0.0 


higher central temperature 



erg/s), but still much less than the total integrated 
energy generation rate in the simulated volume. 
The excess energy deposited on the grid gradually 
leads to an increase in central temperature as de- 
scribed in Section 2. Our anelastic model, how- 
ever, is constructed around a time-independent 
reference state. The inability to follow the evo- 
lution of the background state forces us to exam- 
ine several "snapshots" in the evolution. In this 
paper we consider two different central tempera- 
tures: T Cj 8 = 7.0 and 7.5. Our goal is to estab- 
lish a pseudo-steady-state at each of these epochs, 
which, while being physically unrealistic, might 
still allow us to infer something about the spec- 
trum of temperature fluctuations and the global 
flow pattern. To do this we compensate for the 
excess energy by adjusting the outer luminosity 
to balance the total energy generation rate in our 
computational domain. In order to establish a 
sufficiently large negative entropy gradient, 5s is 
forced to become negative at the outer boundary. 
This is artificial - in a real white dwarf the entropy 
gradient would remain adiabatic throughout the 
convection zone. The compromise was necessary 
in order to achieve a model in near steady state 
that would run stably for a long time. 

In order to allow the model to find a stable so- 
lution quickly, we begin with large turbulent dif- 
fusivities and gradually turn them down, in small 
enough increments that the code is able to find 
a new stable solution. This decrease raises our 
operational Rayleigh and Reynolds numbers. To 
ensure numerical stability the diffusivities must re- 
main large enough to prevent a build-up of entropy 
at the smallest scales. We continue to lower the 
diffusivities to the smallest value that can be ac- 
commodated with a given spatial resolution. We 
then continue to evolve the model at fixed diffusiv- 
ity. In this manner we found pscudo-stcady-states 
for both temperature snapshots. The results pre- 
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Fig. 2. — The horizontally and time averaged ra- 
dial, angular, and total velocities as a function of 
radius in the C7 simulation. The solid line shows 
the root mean square, the dashed lines the mini- 
mum and maximum, and the dotted line (only in 
the total velocity plot) the mean velocities in each 
radial shell. The radial velocity is forced to zero 
at the inner boundary (see text for discussion). 

sented here are from only the last 84, 000 time 
steps, whereas the relaxation phase of the simu- 
lation took about 175, 000 time steps. The long 
duration of the relaxation ensures that any "mem- 
ory" of the initial conditions of the simulation is 
completely erased. 

We considered three models: C7, a non-rotating 
model with a central temperature of T c-8 =7.0, 
C7rot, identical to C7, except that the background 
state is rotated rigidly with an angular velocity 
fl = 0.167 rad/s (Ekman number - 7 x 10~ 4 , 
see Section 4.5), and C75, non-rotating, but with 
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Fig. 3.— Horizontally averaged mass-weighted radial profiles of 5X/X (solid, left ordinate) and 5X (dotted, 
right ordinate), where SX is the entropy, pressure, density, and temperature perturbation (clockwise from 
top left). 



T Ci 8=7.5. The resolution was the same in all three 
models: n max = 237 Chcbyshcv polynomials in 
the radial direction, and spherical harmonics up 
to m max = 42 and ^ max = 85. As we show be- 
low (Section 4), typical plume velocities are ~ 50 
km/s and <~ 100 km/s in the C7 and C75 simula- 
tions, respectively, which translates to convective 

turnover times Of t to = 2(i^ ut er--Rinncr)/«plumc = 

18 and 9 seconds respectively. The two simula- 
tions are evolved for about 70 and 40 seconds star 
time, so for roughly four convective turnovers in 
both cases. The parameters of each simulation are 
summarized in Table 1. 



4. RESULTS 

4.1. Velocities and Thermodynamic Per- 
turbations 

In order for the anelastic approximation to be 
valid, the fluid flow must remain subsonic and 
variations from the reference state must be small, 
<^ 1%. The time and angle averaged (over the 
last 84000 time steps) total, radial, and angular 
velocities are plotted as a function of radius in Fig- 
ure 2. The horizontal mean as a function of radius 
of the absolute value of total velocity increases 
from a minimum of about 50 km/s at R = 400 
km up to about 90 km/s near the center. The 
maximum velocity achieved anywhere on the grid 
is 150 km/s. This is much less than the sound 
speed of c soun( j ^8000 km/s, so the condition of 
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subsonic flow is satisfied. 

Figure 3 shows horizontally averaged mass- 
weighted radial profiles of the thermodynamic per- 
turbations at the end of the C7 simulation. The 
top two plots show the two independent variables, 
the entropy and pressure perturbations. The two 
perturbations are comparable in magnitude, and 
both lie below the 1% level. The maximum en- 
tropy and pressure perturbations in the whole 
volume, of course, are higher: (|<fe|/s) max = 0.015 
and {\SP\/P) max = 0.003. The largest entropy 
perturbation is negative, occurs at the outer 
boundary, and is caused by the large negative 
entropy gradient required to balance the total en- 
ergy generation, as explained in Section 3.3. Using 
Eq. (1) we can determine temperature and den- 
sity fluctuations from the entropy and pressure 
perturbations. These are shown in the lower two 
plots of Figure 3. Since (Hf) p » (MU) and 

/ 01np \ „ / dlnp A 
\dlns J p ^ \dlnPj g > 

^"(M)^''-- (5! " 




1 10 100 

I 



Fig. 4. — Variance of entropy perturbation versus 
spherical harmonic wavenumber I at the end of the 
C7 simulation. The dotted line is a l~ 5 / 3 power 
law, as expected for a Kolmogorov turbulent spec- 
trum. 



For the conditions inside the simulated region, 

(M) « 35 and (fe) s « °- 75 - The lar e est 

positive temperature and density perturbation is 
1.5% and 0.7%, respectively. Since all thermody- 
namic perturbations remain below the few percent 
level, we conclude that the anelastic approxima- 
tion is indeed a valid one for the problem at hand. 

Our resolution of n max = 237, m max = 42, 
and l max — 85 allowed us to lower the diffusiv- 
ities to R — v = 5 x 10 11 cm 2 /s, corresponding 
to Rayleigh and Reynolds numbers of Ra =~ 10 7 
and Re =~ 1500 1 , in all three simulations. The 
transition from laminar to turbulent flow gener- 
ally occurs close to Re=2000, so our simulations 
have not quite reached the regime of fully turbu- 
lent convection, but aren't fully laminar cither. In 
fact, we see about 2/3 of an order of magnitude 
of a Kolmogorov turbulent cascade in an angular 
power spectrum of entropy perturbation, Figure 4. 
The turbulent cascade begins at I = 10 and ex- 
tends down to I = 50. 

While the requirements of a time-independent 
background state and small perturbations prevent 
us from following the increase of the central tem- 
perature all the way to the point where localized 
flames develop, we can look at the C75 snapshot 
to get an idea of what the fluid flow looks like 
at a higher central temperature of 7.5 x 10 8 K. 
This study employed the same resolution as in the 
C7 simulation, but with the more energetic back- 
ground state, we were only able to lower the dif- 
fusivities to R — v = 8.5 x 10 11 cm 2 /s. Due to the 
strong temperature dependence of the nuclear en- 
ergy generation rate, the convection is much more 
vigorous in the C75 simulation. Both the veloci- 
ties and thermodynamic perturbations are larger 
than in the C7 simulation. The mean velocity now 
ranges from 70 to 130 km/s, and the maximum 
velocity on the grid is 276 km/s. These num- 
bers and their temperature scaling agree well with 
the analytic expectations (Woosley, Wunsch, & 
Kuhlcn 2003). Typical thermodynamic perturba- 
tions have also increased: the shell averaged tem- 
perature fluctuation reaches 2% at the innermost 
radial bin, with the largest positive perturbation 
reaching 4.4%. Higher temperature perturbations 
and larger radial velocities increase the likelihood 

We used -u ma x in the calculation of Re. Using u mca n would 
lower this estimate by about a factor of three. 
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of an off-center supernova ignition. This possibil- 
ity is further discussed in Section 4.4. 

4.2. Flow Patterns 

The most striking feature of our calculation is a 
large coherent dipolar circulation. Figure 5 shows 
a three-dimensional representation of this flow - 
two iso-radial- velocity surfaces at v r — ±40 km/s. 
The outflowing surface (blue) is predominantly 
located in one hemisphere, the inflowing surface 
(red), in the other. The dipole nature of the flow 
can also be seen in Figure 6, which shows two- 
dimensional equatorial and meridional slices of the 
temperature fluctuation and the radial component 
of velocity. 

In the non-rotating models there is no preferred 
axis in the calculation. The expansion in spheri- 
cal harmonics might be expected to produce an 
artificial preferred axis along 6 — 0, but the ob- 
served dipole is not aligned with this axis. Neither 
is it aligned with any perturbation in the initial 
conditions. The only explanation is some form of 
"spontaneous symmetry-breaking", growing from 
numerical noise in the calculation. 

Admittedly, we cannot address precisely the na- 
ture of the flow at the center with our method. 
As explained in Section 3.3, we solve the anelas- 
tic hydrodynamics equations subject to a solid in- 
ner boundary condition, which we have placed at 
R = 50 km. The sharp drop in the radial velocity 
component at R < 100 km/s is artificial, caused 
by the solid inner boundary condition. The ac- 
tual radial velocity component at the center can 
be obtained by extrapolation of its value at 100 
km to the origin, as suggested by the large r.m.s. 
angular velocities near the center. This implies 
that the center of the star is not at a calm region 
at all, as one might expect from one-dimensional 
calculations or multidimensional calculations that 
do not carry a full 360 degrees. Similar conclu- 
sion were reached by a recent 2D anelastic study 
of convection in the interior of giant gas planets 
(Evonuk & Glatzmaier 2005). 

4.3. The Temperature Fluctuation Spec- 
trum 

As discussed in Woosley, Wunsch, & Kuhlen 
(2003), the nature of the probability density func- 
tion (PDF) of temperature fluctuations is impor- 



tant for determining whether the supernova ex- 
plosion will have one or multiple ignition points. 
The sharper the high temperature tail of the PDF, 
the more material there is with temperature just 
a bit less than the ignition temperature of the first 
flame. 

We have numerically estimated the tempera- 
ture fluctuation PDF in our simulation for con- 
stant radius shells of thickness 5R ~ 2 km. In each 
shell we subtract from every temperature fluctu- 
ation the mean (ST), divide by the r.m.s. 8T lms , 
and use the resulting quantity as our independent 
variable: (ST - (8T))/ST Tms . We then determine 
the PDF by calculating the fraction of volume of 
the shell occupied by a given fluctuation and divid- 
ing by the total volume of the shell. The resulting 
PDF at three radii (485, 279, and 75 km), cor- 
responding to 20 km from the top, the center of 
the modeled region, and 20 km from the bottom, 
are shown as histograms in Figure 7. We have 
found best-fitting Gaussian distribution (GPDF) 
and exponential distribution (EPDF) and plotted 
them as thin solid and dashed lines, respectively. 
The PDF were fit via simple x 2 minimization, con- 
strained to be centered at ST = {ST), and only the 
positive fluctuations were used in the fit. 

In all three locations the GPDF produced a bet- 
ter fit than the EPDF. None of the fits match 
the distribution perfectly over the entire range 
of positive fluctuations, but in particular in the 
R = 75km shell the Gaussian fit follows the dis- 
tribution out to 3.5<ST rms . The measured distribu- 
tions at R — 485 and 279 km exhibit a strong ex- 
cess of negative temperature fluctuations over ei- 
ther GPDF or EPDF. This excess is entirely com- 
posed of downwards flowing material, making up 
the return current of the global dipole flow de- 
scribed in Section 4.2. 

Niemela et al. (2000) experimentally deter- 
mined the temperature fluctuation PDF in an 
incompressible fluid (cryogenic helium at 4 K) as 
a function of Raylcigh number. At Raylcigh num- 
bers comparable to the ones we have simulated in 
this study, they also found a Gaussian PDF. At 
higher Rayleigh number, however, they observed 
a transition to an Exponential PDF, which re- 
mained intact all the way up to Ra ~ 10 15 , the 
limit of their experiment. Even this is many or- 
ders of magnitude below the Ra ~ 10 25 expected 
in the centers of white dwarfs, and for all practical 
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Fig. 5. — Two iso- velocity surfaces at |u r | = 40 km/s. The sharp division between infalling (red) and 
outflowing (blue) material clearly demonstrates the dipole nature of the flow. 



purposes the true nature of the temperature fluc- 
tuation PDF before the onset of explosive carbon 
burning must be considered unknown. Further 
numerical studies at higher Rayleigh numbers are 
necessary, and may well detect a transition to an 
exponential PDF. 

4.4. The Size and Persistence of Temper- 
ature Fluctuations 

While the PDF captures how rare a tempera- 
ture fluctuation is, averaged over space and time, 
it contains no information about the size and du- 
ration of these fluctuations. The persistence of a 
fluctuation being carried outward is of crucial im- 
portance in determining how far from the center 



the hot spot will run away and ignite the super- 
nova explosion. A hot spot of a given size will be 
shredded to pieces on a timescale set by the degree 
of turbulence within the outflow. Together with 
the laminar outflow velocity, this time scale sets 
the maximum distance from the center that a hot 
spot can reach. As we discussed in Section 4.1, we 
see in our simulations the beginning of a turbulent 
cascade. Unfortunately, our resolution is nowhere 
near adequate to follow this cascade from the inte- 
gral length scale (the size of the largest eddies) all 
the way down to the Kolmogorov scale (the small- 
est eddies). We can, however, make use of scaling 
relations to infer the minimum size fluctuations 
that will survive out to a given distance from the 
center. 
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Fig. 6. — Slices of ST (top) and v r (bottom), in meridional (left) and equatorial (right) planes in the non- 
rotating T Ci 8 = 7 simulation. Solid contours denote positive, dashed ones negative values. 



As discussed by Woosley, Wunsch, & Kuhlen 
(2003), the spots within the dipole flow that run- 
away first will be those in which heating by burn- 
ing just compensates cooling by adiabatic expan- 
sion. This consideration, plus a dipole flow speed 
of approximately 50 - 100 km s" 1 , naturally gives 
ignition about 100 km out from the center of the 
star. 

Within the dipole flow we find a residual ran- 
dom velocity field. The smallest resolved ed- 
dies have a length scale of about 50 km and a 
typical velocity of about 10 km s _1 . Assuming 
Kolmogorov-Obukhov scaling (v oc L 1 / 3 ) we can 



extrapolate to below our resolution limit. We es- 
timate that an eddy 1 km in size will turn over 
in about 1 second. Smaller eddies will turn over 
faster, so we take this as a characteristic volume 
for the fluid that will mix during the 1 to 2 seconds 
it takes a hot perturbation to move to the ignition 
radius - 100 km. 

All these numbers will require a much better re- 
solved calculation to be taken precisely, but they 
suggest that the perturbations that ignite the run- 
away will be moderately large, ~ 1 km, a value 
that might even be resolved in future numerical 
studies. 



11 




-4 -3 -2 -10 1 2 3 4 -3 -2 -10 1 2 3 4 -3 -2 -10 1 2 3 4 
(<5T - <<5T>) / 5T rms (<5T - <<5T>) / <5T rms (<5T - <<5T>) / <5T rms 



Fig. 7. — The probability distribution of temperature fluctuation evaluated in three constant radius shells 
(485, 279, and 75 km), when the central temperature is T c .% = 7.0. The y-axis shows the fraction of the area 
in this shell that intersects fluid elements with temperatures deviating by the amount shown from the mean, 
divided by the root-mean-square of the temperature fluctuations in that shell. The thin solid and dashed 
lines show best- fit Gaussian and Exponential PDF, respectively. Only positive fluctuations were used in the 
fit, and the PDF were constrained to be centered at zero. 



Our results also show significant correlation be- 
tween hot spots. They are not scattered randomly 
throughout the dipole flow, but concentrated near 
its axis. Temperature differences will be ampli- 
fied in the subsequent runaway by the high power 
of the reaction rate temperature sensitivity, but a 
reasonable picture of the ignition might be a sin- 
gle, broad, off-center "plume" rather that multiple 
hot spots. 

4.5. The Effect Of Rotation 

In addition to the non-rotating models dis- 
cussed in the previous sections, we also calculated 
one model (C7rot) with a rigidly rotating back- 
ground state. Rotation is introduced into the sim- 
ulation by solving the anelastic conservation equa- 
tions in a rotating reference frame. This results in 
the Coriolis force term on the right hand side of 
Eq.(2b) of 2p v x f2, where fj is the reference 
state rotation rate vector in rad/sec. The cen- 
trifugal term (p x ft x r) is neglected, because 
it is usually much smaller than the gravitational 
force. 2 

Accreting white dwarfs are thought to rotate 
quite rapidly, with rotation rates approaching Ke- 

2 The centrifugal term must be included for situations where 
the centrifugal and gravitational forces are comparable. 



0.20 




Fig. 8. — Distribution (volume fraction) of the 
ratio of the buoyancy force (— g Sp/p) to the ra- 
dial component of the Coriolis force (2 (v x fl) r ). 
Material to the left of the dotted vertical line ex- 
periences a Coriolis force that is stronger than the 
local buoyancy force. 
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Fig. 9. — Like Figure 6, but for the C7rot simulation, 
dipole flow. 

plerian at the surface (eg. Yoon & Langer 2004). 
In the KEPLER model which is the basis for this 
anelastic study, the Keplerian rotation rate at the 
surface is Cl = {GM/R 3 ) 1 / 2 = 5.6 rad/sec. We 
are interested in determining what effect even a 
small amount of rotation has on the dipole flow 
that we observed in the non-rotating simulations. 
For this purpose, we chose a moderate background 
rotation rate of 0.167 rad/sec, only 3% of fio- As 
we shall see, even this comparatively small amount 
of rotation has an appreciable effect upon the fluid 
flow. At the outer boundary of our simulation 
the ratio of centrifugal to gravitational force is 




in:i 



A rotating background state seems to break up the 



n 2 R 3 /(GM) = 7 x 10~ 5 , justifying the omission 
of the centrifugal force term. The importance of 
the coriolis force relative to buoyancy forces is ad- 
dressed in Figure 8. We have plotted the fraction 
of the total simulated volume at a given ratio of 
buoyancy force (—gSp/p) to the radial component 
of coriolis force versus the logarithm of this ratio. 
This plot shows that buoyancy dominates corio- 
lis forces for most of the simulated volume: only 
for 1.3% of the total volume does the radial com- 
ponent of the coriolis force exceed the buoyancy 
force. Two further dimensionless numbers char- 
acterize the importance of rotation. The Ekman 
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number is a measure of the relative importance 
of viscous to Coriolis, and the Rossby number of 
inertial to rotational forces. 



Ek 



2QD 2 



Ro 



(6a) 



(6b) 



2Dil sin 9 ' 

where 9 is the co-latitude. We have Ek w 7 x 1CT 4 , 
and at the equator Ro=l or 1/3, depending on 
whether we use w max or w me an- 

Even though the centrifugal and coriolis forces 
are small compared with gravity and buoyancy, 
respectively, the rotating background state signif- 
icantly alters the flow in the simulation. Figure 9 
shows equatorial and meridional slices of tempera- 
ture perturbation and radial velocity for the C7rot 
simulation. A comparison with Figure 6 shows 
that the dipole flow has pretty much disappeared. 
Instead of a well defined outflow (inflow) in the 
lower right (upper left) quadrant of the meridional 
slice in the C7 simulation, both inflowing and out- 
flowing plumes can be found in all four quadrants 
of the same slice in the C7rot simulation. The tem- 
perature fluctuations in C7rot are also broken up 
into multiple hot and cold spots, instead of one hot 
and one cold plume in C7. Two-dimensional simu- 
lations at a much higher Raylcigh number show a 
transition from a dipole to a differentially rotating 
longitudinal flow as the rotation rate is increased 
(Evonuk & Glatzmaier 2005). It is important to 
note, however, that flow through the center might 
still be possible in a different configuration. The 
temperature fluctuations and convective velocities 
are of the same order as in the C7 simulation, and 
the results of Section 4.3 and 4.4 will continue to 
hold as long as there is at least one central outflow 
of hot material. An asymmetric supernova explo- 
sion, however, becomes less likely in the absence 
of a strong dipole flow. 

5. CONCLUSIONS 

Three-dimensional studies at moderate Rayleigh 
number (~ 10 7 ) of carbon ignition in the highly 
degenerate core of a white dwarf star show a dipole 
nature to the flow. If such a flow pattern persists 
at the much higher Rayleigh number of the actual 
white dwarf, this may lead to the asymmetric ig- 
nition and explosion of Type la supernovae. We 



also find that a moderate degree of rotation, corre- 
sponding to less than 10% critical at the surface, 
disrupts this dipole, leading to ignition over a 
broader range of angles. 

Within the flow, there still exists turbulence 
and this will limit the persistence of temperature 
fluctuations that might serve as ignition points. 
Only the larger fluctuations will survive the tran- 
sit through the core to the estimated ignition ra- 
dius of 100 km. We estimate those fluctuations 
will be a km or larger. 

Finally we have estimated the PDF for the tem- 
perature fluctuations (Gaussian) and commented 
on the degree of spatial correlation for the hot 
spots (high). This suggests that once a runaway 
ignites in some locality, other points will swiftly ig- 
nite in approximately the same vicinity (although 
the ignition may continue for some time). 

Though the present study has been the first to 
follow the burning in 3D in a large fraction of the 
unstable core, it has raised many questions that 
will require further study. 

Chief among them are the scaling properties 
of the dipole flow with Reynolds number. The 
present study has a Reynolds number approxi- 
mately 10 11 times smaller than the actual star. It 
is quite likely that the actual flow becomes increas- 
ingly chaotic at larger Reynolds number, analo- 
gous to what is seen in Rayleigh-Bcrnard convec- 
tion (Kadanoff 2001). Two-dimensional simula- 
tions are capable of showing the dipole flow when 
run for full cylinders, and might have adequate 
resolution to study this scaling. 

Second, though we have argued that its effects 
are minimal, the artificial inner boundary condi- 
tion in the present study is troublesome and needs 
to be removed. This is not an inherent deficiency 
in the anelastic method, but results from the use of 
a spectral method in spherical coordinates - the 
angular derivatives diverge at r = 0. A Carte- 
sian, finite-volume version of the code has been 
developed by Evonuk & Glatzmaier (2005) and is 
currently being used to further study the present 
problem. Initial 2D simulations of convection in 
giant gaseous planets with this code show that the 
presence of a solid core in the non-rotating case 
can lead to a significantly different flow structure. 
With increasing rotation rate, however, the flow 
settles into a differentially rotating profile that is 
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only weakly affected by the presence of a solid 
core. Thus it may not be necessary to remove the 
solid inner core for simulations of rapidly rotating 
white dwarfs. 

Third, the study needs to be started earlier (at 
a lower central temperature) and run longer (very 
many convective turnover times). It takes time for 
the ID background state model to adjust to the 
new code. Ideally, one wants to watch the gradual 
rise in the overall temperature until the runaway 
actually ignites in localized regions. One can do 
that by gradually adjusting the background state 
in the anelastic model. Here we sought stability at 
the expense of imposing an artificially large radia- 
tive boundary at the outer edge of the problem. 

Given the importance of the "ignition problem" 
to understanding Type la supernovae, we expect 
significant progress on all these fronts in the near 
future. 

We thank M. Zingale for informative discus- 
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NASA (NAG5-12036), and the DOE Program for 
Scientific Discovery through Advanced Comput- 
ing (SciDAC; DE-FC02-01ER41176. All computa- 
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processor IBM SP RS/6000 supercomputer. 

REFERENCES 

Baraffe, I., Hcger, A., & Woosley, S. E. 2004, ApJ, 
615, 378 

Boyd, J. P. 2000, Chebyshev and Fourier Spectral 
Methods (2nd ed.; Mineola, NY: Dover) 

Evonuk, M., & Glatzmaier, G. A. 2005, Planetary 
and Space Science, submitted 

Gamezo, V. N., Khokhlov, A. M., Oran, E. S., 
Chtchelkanova, A. Y., & Rosenberg, R. O. 2003, 
Science, 299, 77 

Garcia-Senz, D., & Woosley, S. E. 1995, ApJ, 454, 
895 

Garcfa-Senz, D., & Bravo, E. 2005, A&A, 430, 585 

Glatzmaier, G. A. 1984, J. Comp. Phys., 55, 461 

Glatzmaier, G. A., & Roberts, P. H. 1995, Nature, 
377, 203 



Hoflich, P., & Stein, J. 2002, ApJ, 568 779 

Hillebrandt, W., & Nicmeycr J. 2000, ARAA, 38, 
191b 

Kadanoff, L. P. 2001, Physics Today, 54(8), 34 

Kraichnan, R. H. 1962, Phys. of Fluids, 5, 1374 

Lantz, S. R., & Fan, Y. 1999, ApJ, 121, 247. 

Niemela, J. J., Skrbek, L., Sreenivasan, K. R., & 
Donnelly, R. J. 2000, Nature, 404, 837 

Niemeyer, J. C, & Hillebrandt, W. 1995, ApJ, 
452, 769 

Niemeyer, J. C, Hillebrandt, W., & Woosley, S. 
E. 1996, ApJ, 471, 903 

Niemeyer, J. C, & Woosley, S. E. 1997, ApJ, 475, 
740 

Niemeyer, J. C. 1999, ApJL, 523, L57 

Hillebrandt, W., & Niemeyer, J. C. 2000, 
ARA&A, 38, 191 

Nomoto, K., Sugimoto, D., & Neo, S. 1976, 
Ap&SS, 39L, 37 

Nomoto, K., Thielemann, F.-K., & Yokoi, K. 1984, 
ApJ, 286, 644 

Plewa, T., Calder, A. C, & Lamb, D. Q. 2004, 
ApJ, 612, L37 

Porter, D. H., & Woodward, P. R. 2000, ApJS, 
127, 159 

Ropke, F. K., & Hillebrandt, W. 2005, A&A, 431, 
635 

Ropke, F. K., Hillebrandt, W., Niemeyer, J. C, & 
Woosley, S. E. 2005, A&A, submitted 

Sun, Z.-P., Schubert, G., & Glatzmaier, G. A. 
1993, Science, 260, 661 

Timmcs, F.X., & Swesty, F.D., 2000, ApJ Supple- 
ment Series, 126, 501 

Weaver, T. A., Woosley, S. E., & Zimmerman, G. 
B. 1978, ApJ, 225, 1021 

Woodward, P. R., Porter, D. H., & Jacobs, M. 
2003, ASP Conf. Ser. 293: 3D Stellar Evolution, 
293, 45. see also www.lcse.umn.edu/3Dstars. 



15 



Woosley, S. E., Wunsch, S., & Kuhlcn, M. 2003, 
ApJ, 607, 921 

Wunsch, S., & Woosley, S. E. 2004, ApJ, 616, 1102 

Yoon, S.-C, & Langer, N. 2004, A&A, 419, 623 



This 2-column preprint was prepared with the AAS IATfrjX 
macros v5.2. 



16 



